
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Previously saved workspace restored]

> rm(list=ls())
> gc()
         used (Mb) gc trigger (Mb) max used (Mb)
Ncells 331916 17.8     592000 31.7   424061 22.7
Vcells 420573  3.3    1023718  7.9  1012053  7.8
> require(foreign)
Loading required package: foreign
> require(ineq)
Loading required package: ineq
Warning message:
package 'ineq' was built under R version 3.3.2 
> require(dplyr)
Loading required package: dplyr

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

> library(Hmisc)
Loading required package: lattice
Loading required package: survival
Loading required package: Formula
Loading required package: ggplot2

Attaching package: 'Hmisc'

The following objects are masked from 'package:dplyr':

    combine, src, summarize

The following objects are masked from 'package:base':

    format.pval, round.POSIXt, trunc.POSIXt, units

> library(lavaan)
This is lavaan 0.5-22
lavaan is BETA software! Please report any bugs.
Warning message:
package 'lavaan' was built under R version 3.3.2 
> 
> 
> setwd(paste(substr(getwd(),1,regexpr("Dropbox",getwd())[1]-1),"Dropbox/coauthors/AMT_Bisbee_Larson/bisbee_larson_replication_final",sep=""))
> getwd()
[1] "D:/Dropbox/coauthors/AMT_Bisbee_Larson/bisbee_larson_replication_final"
> source("./tools/leg2.R")
> ###################################################################################
> ##
> ##  File Name:    protest_results.R
> ##
> ##  Input Files:  prot-simple.dta
> ##  Output Files: SI_prot_causal.pdf             (Figure 5)
> ##
> ##  Purpose:  This R script summarizes the causal relationship between tie dimension
> ##            and decision to join a protest by exploiting random assignment of 
> ##            elicited tie in the data. 
> ##
> ##############################################################################################
> 
> ttest.coefs <- function(beta1,beta2,se1,se2) {
+   t.stat <- (beta1-beta2)/(sqrt(se1^2 + se2^2))
+   pval <- 2*pt(-abs(t.stat),df = 450)
+   star <- paste("(",sprintf("%.2f",round(pval,2)),ifelse(pval < 0.05,"*)",")"),sep="")
+   return(star)
+ }
> 
> 
> # ----------------------------------------- Causal ------------------------------------- #
> dat.cause <- read.dta("./protests/prot-simple.dta")
> dat.cause$prot_dim_interaction <- dat.cause$prot_dim_interaction*-1
> dat.cause$prot_dim_seat_synth2_mst <- dat.cause$prot_dim_seat_synth2_mst*-1
> dat.cause$prot_dim_seat_synth2_mst_rob <- dat.cause$prot_dim_seat_synth2_mst_rob*-1
> 
> dims <- c("jobsearch","contribute","garner","personalgain","grpgain","homo_pol","homo_rel","homo_educ",
+           "homo_ls","pers_cris_synth","pers_succ_synth","prof_cris_synth","prof_succ_synth",
+           "seat_synth2_lst","seat_synth2_mst","interaction")
> coefs.biv <- ses.biv <- NA
> for(i in 1:length(dims)) {
+   coefs.biv[i] <- summary(lm(paste("prot_join ~ prot_dim_",dims[i]," + prot_peace",sep=""),dat.cause))$coefficients[2,1]
+   ses.biv[i] <- summary(lm(paste("prot_join ~ prot_dim_",dims[i]," + prot_peace",sep=""),dat.cause))$coefficients[2,2]
+ }
> names(coefs.biv) <- names(ses.biv) <- dims
> 
> model.cause <- paste("prot_join ~ ",paste(paste("prot_dim_",dims,sep=""),collapse=" + ")," + weakest_dum + weak_dum + strong_dum + strongest_dum + prot_peace",sep = "")
> fit.cause <- sem(model.cause,data = dat.cause)
> coefs.cause <- fit.cause@ParTable$est[which(grepl("^~$",fit.cause@ParTable$op))][1:16]
> ses.cause <- fit.cause@ParTable$se[which(grepl("^~$",fit.cause@ParTable$op))][1:16]
> 
> ords <- order(coefs.biv)
> labs <- c("Job Search","Contrib. to Entrep.","Garner Conts.","Personal Gain (% $100)","Group Gain (% $100)","Political Hom.",
+           "Religious Hom.","Education Hom.","Class Hom.","Pers. Crisis","Pers. Success",
+           "Prof. Crisis","Prof. Success","Least in Common","Most in Common","Pref. Interaction")
> 
> pdf("./1_Figures/SI_prot_causal.pdf")
>   nf <- layout(matrix(c(1,2,3),1,3,byrow=T),widths=c(2,3,3))
>   #layout.show(nf)
>   par(mar=c(2,5,3,0))
>   plot(rep(0,length(coefs.cause)), 1:length(coefs.cause), 
+        type = "n", axes = F, xlab = "", ylab = "",main = "")
>   axis(2, at = 1:length(coefs.cause), label = labs[ords], las = 1, tick = FALSE, cex.axis =.9,pos = .9)
>   
>   par(mar=c(2,0,3,1))
>   pchs <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,19,21)
>   cexs <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,.6,1.2)
>   cols <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,rgb(0,0,0,.1),rgb(0,0,0,.8))
>   ltys <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,1,1)
>   bgs <- ifelse(abs(coefs.biv[ords]/ses.biv[ords]) < 2.33,NA,"white")
>   plot(coefs.biv,1:length(coefs.biv),type='n',xlab="Coefficient",ylab="",yaxt='n',xlim=c(-.1,.2),main="Bivariate Regressions")
>   segments(coefs.biv[ords] - 2*ses.biv[ords],1:length(coefs.biv),coefs.biv[ords] + 2*ses.biv[ords],col=cols,lty=ltys)
>   points(coefs.biv[ords],1:length(coefs.biv),pch=pchs,col=cols,cex=cexs,bg=bgs)
>   abline(v=0,lty=2)
>   
>   pchs <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,19,21)
>   cexs <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,.6,1.2)
>   cols <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,rgb(0,0,0,.1),rgb(0,0,0,.8))
>   ltys <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,1,1)
>   bgs <- ifelse(abs(coefs.cause[ords]/ses.cause[ords]) < 2.33,NA,"white")
>   plot(coefs.cause,1:length(coefs.cause),type='n',xlab="Coefficient",ylab="",yaxt='n',xlim=c(-.15,.15),main="SEM Analysis")
>   segments(coefs.cause[ords] - 2*ses.cause[ords],1:length(coefs.cause),coefs.cause[ords] + 2*ses.cause[ords],col=cols,lty=ltys)
>   points(coefs.cause[ords],1:length(coefs.cause),pch=pchs,col=cols,cex=cexs,bg=bgs)
>   abline(v=0,lty=2)
> dev.off()
null device 
          1 
> 
> proc.time()
   user  system elapsed 
   1.79    0.35    2.14 
